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1.  Introduction 


The  problem  of  axisymmetric  indentation  of  an  elastic  half-space  by  a  rigid  flat-ended  cylindrical 
punch  (the  Boussinesq  problem)  is  summarized  in  this  report.  Exact  analytical  solutions  to 
elastic  indentation  problems  are  important  for  verification  of  numerical  codes  capable  of 
modeling  such  phenomena  (2,  3),  and  for  modeling  the  onset  of  permanent  deformation  and 
fracture  (5),  a  topic  of  active  research  in  the  Materials  and  Manufacturing  Sciences  Division  of 
the  U.S.  Army  Research  Laboratory. 

Boussinesq  (6)  and  Love  (7)  were  among  the  first  to  consider  such  a  problem,  but  Sneddon  (4,  8) 
solved  the  mixed  boundary  value  problem  for  the  rigid  flat-ended  punch  in  cylindrical  polar 
coordinates  by  using  the  Hankel  transform  to  convert  the  governing  equation,  i.e.,  the  biharmonic 
equation  V4^  =  0,  into  a  set  of  dual  integral  equations  solvable  by  methods  introduced  by 
Titchmarsh  (9).  The  generality  of  Titchmarsh’s  methods  (9)  were  such  that  they  were  also  used 
by  Busbridge  (10),  and  Tranter  (11)  to  derive  solutions  to  the  electrified  disk  problem  important 
in  potential  field  theory,  V2©  =  0.  Other  powerful  methods  of  analysis  rely  on  using  the 
Papkovich-Neuber  solution  (Galerkin  vector),  which  for  body-force  free  axisymmetric  problems 
involves  determining  two  scalar  harmonic  functions  solvable  by  Hankel  transforms  (12).  Hankel 
transforms  and  dual  integral  equations  have  been  used  extensively  by  a  number  of  researchers  to 
solve  axisymmetric  boundary  value  problems  in  fluid  (13),  poroelastic  (14),  piezoelectric  (15), 
and  magnetoelectroelastic  (16)  media. 

We  focus  our  study  on  the  fully  three-dimensional,  axisymmetric  indentation  problem  despite  the 
long  history  of  developments  in  the  two-dimensional  theory  of  mixed  boundary  value  problems, 
known  as  Riemann-Hilbert  problems,  and  their  solution  via  the  Plemelj  (17)  formulae;1 
two-dimensional  half-plane  solutions  for  rigid  flat  punch  indentation  problems  are  decidedly 
non-physical,2  since  vertical  surface  displacements  become  unbounded  away  from  the  punch. 

This  result  has  been  obtained  by  both  England  (27;  page  70,  equation  (3.67))  using  complex 
variable  methods,  and  Sneddon  (7;  page  433,  equation  (115))  using  Lourier  transform  methods. 

During  the  course  of  our  literature  survey  on  axisymmetric  indentation,  we  encountered  many 
typographical  errors  or  misprints  of  equations  that  have  been  unwittingly  propagated  through  the 

barton  and  Perlin  (18)  refer  to  these  equations  as  the  Sokhotskii-Plemelj  formulae,  and  Gakhov  (19)  refers  to 
these  equations  as  the  Sokhotskii  formulae,  since  Sokhotskii  derived  them  first  in  his  1873  dissertation. 

Solutions  to  boundary  value  problems  that  are  non-physical  may  be  more  useful  for  code  verification  than  exact 
physical  solutions  with  field  singularities;  see,  e.g.,  the  method  of  manufactured  solutions  in  Oberkampf  and  Roy  (20). 
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literature  by  various  authors;  we  have  identified  these  errors  in  this  report  for  readers  interested  in 
such  analysis.  Furthermore,  the  original  works  of  Sneddon  (1,  4,  8)  lacked  analytical  expressions 
for  the  vertical  displacements  within  the  half-space  region;  indeed,  only  expressions  for  the 
vertical  displacements  of  the  surface  of  the  half-space  are  given,3  and  we  wondered  why  the  more 
general  half-space  solution  was  omitted  by  Sneddon.  The  vertical  displacement  solution  in  the 
entire  half-space  is  important  for  verification  of  numerical  codes.  We  traced  the  probable  cause 
of  the  omission  to  the  unevaluated  J°  term  appearing  in  the  expression  for  the  vertical 
displacement  uz.4  In  general,  the  terms  can  be  defined  (7)  by  taking  the  imaginary  S  part  of 
a  Laplace  transform  involving  mth- order  Bessel  functions  Jm(pp ) 


noo 

Jn=%  Pn-1e-p«-i>Jm(pp)dp  .  (1) 

Jo 

Explicit  solutions  to  equation  1  can  be  obtained  by  taking  the  imaginary  part  of  a  more  general 
result  found  in  Sneddon’s  textbook  (1\  page  514,  appendix  A), 


■C  =  3  L-”p”(C  -  +  n)  2Ft  (^±4  i(m  +  a  +  1);  m  +  1;  -^A_)  )  ,  (2) 

where  ( p ,  Q  >  0, 5R(m  +  n)  >  0,  and  2F\  is  the  regularized  hypergeometric  function 
■2Fi  =  2Fi(a,  b ;  c;  z) /T(c).  According  to  Watson  (22;  page  384),  both  Hankel  (23)  and 
Gegenbauer  (24)  “proved  that  5R(m  +  n)  >  0,  to  secure  convergence  at  the  origin.”  More  direct 
is  the  observation  that  for  the  Jq  term,  if  m  —  n  —  0,  then  r(m  +  n)  =  (m  +  n  —  1)!  is  undefined 
in  equation  2.  We  also  found  in  the  Tables  of  Integral  Transforms  from  the  Bateman  Manuscript 
Project  (25;  page  182,  equation  (5))  that  a  solution  to  the  singular  integral  (equation  1)  does  not 
exist  for  m  —  n  —  0, 


JP  =  9 


e  p(c  Jmjpp )  dp  =  ^t  ^  ^ 

/o  V 


—m  pin 


1  +  A  / 1  + 


(C-d2 


m 


m  >  0 


(3) 


However,  from  the  same  work  (25),  one  can  find  an  explicit  representation  for  the  Jq  term  given 
3See  e.g.,  Sneddon’s  textbook  (7;  page  461,  the  equation  is  not  numbered). 

4See  e.g.,  Sneddon’s  textbook  (7;  page  459,  equation  (44)),  and  reproduced  for  convenience  in  equation  5  of  this 
report. 
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on  page  101,  equation  (17)  in  terms  of  a  Fourier  sine  transform, 


7°  — 
Jo  — 


"°°  e"K 


sin  (p)Jo(pp) 
P 


dp  =  sin 


-l 


\/c2+(p-i)2+vc2+(p+i)2 


(4) 


Thus,  a  formal  relationship  should  be  derivable  between  the  solution  for  Jq  given  in  equation  4 
and  equation  3  evaluated  at  m  =  0;  we  derive  such  a  relationship  in  this  report  through 
application  of  L’Hopital’s  rule  to  equation  3.  In  fact,  we  derive  new  expressions  for  all  the  .1"' 
terms  that  appear  as  stress  and  displacement  components  in  Sneddon’s  solution  to  the 
axisymmetric  indentation  problem  (7).  The  new  analytical  solutions  are  used  to  verify  our 
peridynamics  numerical  code  (2,  3),  and  illustrate  principal  stress  trajectories  in  an  elastic 
half-space  under  the  action  of  a  rigid  cylindrical  punch.  Sneddon’s  classic  textbook  (7)  was 
published  only  three  years  before  the  publication  of  the  Bateman  treatise  (25)  so  the  reason  why 
Sneddon  was  not  aware  of,  or  did  not  publish  the  result  given  by  equation  4,  remains  an  open 
question  and  was  in  fact  another  motivation  for  this  research. 


2.  The  Mechanics  of  Indentation 


The  solutions  for  the  displacements  and  stresses  in  an  elastic  half-space  under  the  action  of  a  rigid 
flat-ended  cylindrical  punch  that  appear  in  Sneddon’s  textbook  (7)  are  listed  below: 


(Tg  = 


Ur  = 


2pe  (Ji  - 

7t(A  +  2  fi) 

2‘  +  Jg 


u7  = 


7T 


CTZ  = 


4/re  ((J®  +  Ji)  (A  +  /i) 
7TCl(  A  +  2  n) 

+  /i) 


na{\  +  2/r) 

4J»AM£  V*  (jo  -  Ji) 


<7  r  -\~  <7  0  — 


7ra(A  +  2/r)  nap(\  +  2/r) 

4/re  (^/i  (2A  +  /r)  —  C^2*(^  +  d)) 


7ro(A  +  2/r) 


(5) 
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In  these  equations,  the  dimensionless  radial  coordinate  p  —  r/a  is  the  physical  radial  distance  r 
normalized  by  the  indenter  (punch)  radius  a,  and  the  dimensionless  depth  coordinate  (  =  z/a  is 
the  physical  depth  z  normalized  by  the  punch  radius  a.  e  is  the  indentation  depth,  and  the  Lame 
parameters  are  given  by  A  and  /j.  The  radial  stress  ay  is  not  given  in  Sneddon’s  textbook  (I)  but 
can  be  derived  readily  by  subtraction  of  ay  from  ay  +  ay  from  equations  5  resulting  in, 


O'  r  — 


Ape(J]p  —  (A  +  /i)(p(J”  —  (J$)  +  C^i1)) 


nap(\  +  2/i) 

where  the  J™  terms  are  defined  as  singular  integrals  involving  Bessel  functions  Jm(pp), 


(6) 


roo 

Jn  =  /  pn~1e^psm(p)Jm(pp)dp  ,  (7) 

Jo 

or  as  Laplace  transforms  involving  Bessel  functions,  i.e.,  the  imaginary  part  of  the  singular 
integral, 


POO 

•C  =  »  /  pn-le-^Jm{pp)dp  .  (8) 

Jo 

Explicit  analytical  expressions  for  a  number  of  specializations  of  the  integral  in  equation  8  can  be 
found  in  the  Tables  of  Integral  Transforms  from  the  Bateman  Manuscript  Project  (25;  page  182; 
equation  (5))  for  example,  if  m  >  0  and  n  =  0, 


J0m  =  3  /  p-le-p(^:Klrn(pp)  dp  =  Z 

Jo 

Thus,  for  m  =  1,  equation  9  becomes 


( (C  -  (i  + 


m 


(9) 


Jo=% 


P 


(C-t)(l  +  H(p,C)) 


(10) 
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where,  H(p ,  ()  =  yl  +  .  However,  for  m  =  0,  care  must  be  taken  for  evaluation  of  the 

integral  in  equation  9  via  a  limit  process,  i.e., 


7°  — 
J  o  ~ 


lim 


jm  _ 

J0  — 


/M 

lim  — — 

gym) 


=  lim  (r 

m—>0 


(C-0 


-m  pm(i+^(p,or 


m 


(ID 


In  evaluation  of  the  limit  in  equation  11,  since  /(m)  =  (v(l)  =  g(m )  =  0  as  5ft(m)  — >  0,  and 
since  the  limit  of  the  ratio  f'(m)/g'(m)  as  'R(rn)  — »  0  exists,  we  may  employ  L’Hopital’s  rule  for 
the  evaluation  of  , 


J?  =  lim  ^ 

m->-o  gym) 


f\rn) 
m-y 0  g'(m ) 


=  lim 


(12) 


The  derivative  of  the  numerator  and  denominator  of  equation  1 1  results  in 


r 


-9f((C  ~i)  mpm  (log(l  +  77(p,C))  +  log(C 

lim  - 


—  l°g(p))  (1  +  H(p,  ())  m) 


•  (13) 


Thus, 


Jo  =  — 0=  (l°g  (T  +  H(p,Q)  +  log(C  ~i)  ~  log(p))  ,  (14) 


or 


7°  —  —A 
j0  —  'J 


log 


(C-»)(i  +  ff(p,0) 

P 


(15) 


Alternatively,  since  z  =  |z|  e*arg^,  then  log(z)  =  In  |z|  +  i  arg(z),  and  9f(log(z))  =  arg(z)5 
resulting  in 

Multivalued  arg(z)  is  defined  by  arg(z)  =  {Arg  (z)  +  2nk  :  kG  Z},  —  tt  <  Arg  (z)  <  tt  . 
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2  =  W  eiarg(z)  =  \z\  eie 


then  equation  1 1  can  be  written  succinctly  as 


lim  S  ( — )  =  lim  9\z\m^^-^-  =  lim  9 1 z\m sine (mO)  =  6  =  arg(^) 
m— >0  V  m  )  to— >0  m6  m— >0 


Interestingly,  we  later  became  aware  of  an  expression  for  Jq  in  Barquins  and  Maugis  (27;  page 
355;  appendix  2),  who  obtained  it  from  an  earlier  work  of  Dahan6  (28),  although  the  steps  taken 
for  its  derivation  do  not  appear  in  either  work.  The  expression  that  appears  in  Barquins  and 
Maugis  (27)  using  Sneddon’s  ( 1 )  notation  is 


Dahan  (28)  refers  to  the  earlier  work  of  Gerrard  and  Harrison7  (30),  where  in  fact,  an  expression 
for  Jq  can  be  extracted  by  specializing  their  general  results  for  an  anisotropic  elastic  half-space  to 
an  isotropic  elastic  half-space  under  vertical  displacement  by  a  rigid,  flat-ended  cylindrical  punch 
(see  e.g.,  pages  12,  14,  and  15  of  Gerrard  and  Harrison  (30)).  The  expression  for  that  we  have 
extracted  from  Gerrard  and  Harrison  (30)  using  Sneddon’s  (1 )  notation  is 

6Dahan  (28)  refers  to  his  thesis  (29)  as  a  possible  source  for  derivation,  which  is  inaccessible  to  us. 

7Gerrard  and  Harrison  (30)  refer  to  an  earlier  reference  from  Koning  (31)  where  we  find  all  of 
Sneddon’s  (1)  expressions  (27)  listed  without  reference;  in  addition,  a  correct  expression  for  Jq  = 

sin^1  ( \l (3  +  p2  +  1  —  \! (C2  +  p2  +  l)2  —  Ap2 / (■s/2p)  \  appears  without  derivation. 
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2 


(20) 


Jq  =  sin 


-l 


Vc + (p  - 1)2  +  vc2  +  ~p  + 1)2 


Equation  20  is  identical  to  that  found  in  the  Bateman  Manuscript  Project  (25),  i.e.,  equation  4 
herein,  and  subsequently  found  in  Fabrikant  (52;  page  221,  equation  (54)).  Fabrikant  (55)  also 
applied  the  theory  of  the  potential  to  solve  a  number  of  mixed  boundary  value  problems  in 
mechanics,  including  the  bonded/unbonded  cylindrical  punch  problem  in  both  isotropic  and 
transversely  isotropic  elastic  media  (52).  Using  the  identity 

sin"1  =  §  —  tan"1(a;),  x  >  0  one  can  verify  that  equation  19  is  identical  to  equation  20. 

Our  own  result  for  Jq,  given  by  equation  15  or  equation  16,  can  be  shown  to  be  equivalent  to  that 
of  Gerrard  and  Harrison  (30)  in  equation  20  using  an  identity  from  Churchill’s  text  (54;  page  62, 
section  29),  that  relates  the  inverse  sine  function  to  a  logarithmic  function, 


sin  1(z)  =  —i  log  ^\/l  —  z 2  +  iz'j 


=  —i  log 


(C  —  *)  (l  +  \J  1  + 


(C-92 


P 


(21) 


and  if  we  solve  equation  21  for  z,  we  find  that  z  = 


_  _i+*C 


.  Thus,  our  result  in  equation  16  follows 


readily  by  simplification  of  sin 


-l 


1+iC 


.  A  plot  of  the  Jq  term  is  illustrated  in  figure  1. 


Barquins  and  Maugis  (27)  have  identified  several  typographical  errors  in  the  literature  by  noting 
“...that  J°  and  ctq  are  misprinted  in...”  Sneddon’s  paper  (4),  i.e.,  equations  (3)  and  (16)  in  that 
work8;  we  also  discovered  another  error  in  the  term  for  ur  appearing  in  equation  (8)  of  Sneddon’s 
paper9  (4);  this  error  was  also  incorrectly  transcribed  to  appendix  2,  equation  (34)  of  Barquins 
and  Maugis  (27).  These  typographical  errors  are  corrected  in  Sneddon’s  textbook  (1). 


Proceeding  in  a  similar  fashion  by  using  the  results  from  (25;  page  182,  equation  (7))  for 
m  —  1,  n  —  2,  equation  8  becomes 

8  Jq  should  be  ./(  in  the  first  term,  and  Jj  should  be  in  the  last  term  of  equation  (3)  in  Sneddon’s  paper  (4);  the 
transcendental  argument  sin  should  appear  as  sin  (4p  —  0)  in  equation  (16)  of  Sneddon’s  paper  (4). 

9The  (Jq  —  (2+U^ji)  term  should  appear  as  (Jq  —  ^ji)  jn  equation  (8)  of  Sneddon’s  paper  (4). 
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Figure  1.  Plot  of  the  Jq  term  given  by  equation  16,  19,  or  20. 


Jl  =  Of 


pe  p ^  l\j\  (pp)  dp  — 


3( 


2pr  (|) 


V"2  <P2  +  K  “  ()2) 


3/2  ■ 


=  9 


p 


((-*)3%()3 


(22) 


Thus,  using  the  result  from  the  Bateman  Manuscript  Project  (25;  page  182,  equation  (1))  for 
m  —  1,  n  —  1,  equation  8  becomes 


Jl  =  %  e-p(C-0  Ji^pp)  dp  = 


P 


((-P(A()(%()  +  1) 


(23) 


and  using  the  result  from  the  Bateman  Manuscript  Project  (25;  page  182,  equation  (1))  for 
m  =  0,  n  —  1,  equation  8  becomes 


^  =  afe^Jo(w)^  =  3(_^__)  .  (24, 

Finally,  using  the  result  from  the  Bateman  Manuscript  Project  (25;  page  182,  equation  (2))  for 
m  —  0,  n  —  2,  equation  8  becomes 
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*  =  *£*-«■**<»>)*>  =  g  ((C  -  .ygfaO’)  '  <M) 

Our  results  for  the  ,7™'  terms  necessary  for  determining  the  stress  and  displacement  fields 
(equations  5  and  6)  in  an  elastic  half-space  under  the  action  of  a  rigid  flat-ended  cylindrical  punch 
are  written  succinctly  as 


The  results  summarized  in  equation  26  are  derived  using  the  Tables  of  Integral  Transforms  from 
the  Bateman  Manuscript  Project  (25).  If  we  use  instead.  Table  VII.  Hankel  Transforms  on  page 


528  of  Sneddon’s  textbook  (/),  we  can  verify  the  results  summarized  in  equation  26.  Using  the 
Hankel  transform  tables  found  in  Sneddon’s  textbook  (7),  however,  we  find  that  the  two  terms  7(j 


and  J\  differ  in  form,  but  are  equivalent  to  those  appearing  in  equation  26.  These  are  given  by 

J.  =  3  ((C-iXyi-D)  =  im  _  miM)  _  1})  and  ji  =  3  (4|!)  =  _I3 


Our  expressions  in  terms  of  dimensionless  coordinates  p  and  Q  differ  in  form,  but  are 
equivalent10  to  those  given  in  Sneddon’s  textbook11  ( 1 )  and  Barquins  and  Maugis  (27),  viz., 

10We  have  verified  this  numerically,  to  machine  precision,  by  comparing  all  terms  in  equation  26  and  equa¬ 
tion  27. 

1  'All  equations  27  appear  originally  in  Sneddon’s  textbook  (7),  except  for  the  Jq  term,  which  appears  explicitly  in 
Barquins  and  Maugis  (27).  The  J\  term  is  misprinted  in  Sneddon’s  textbook  (7;  page  462,  equation  (56)),  as  it  omits 
the  p  that  multiplies  the  sin  term. 
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7° 

Jo 


7T 

- tan 

2 


tO  =  Shl  (f ) 

Vr 

rsm('f-e) 
J'2  ~  RW 

Tl  Tsin(0-f) 

pVR 

psin(f) 
2  i?3/2 

_  1-  y^sin(f) 
°  P 

^C2  +  p2  +  R-l\ 

V2  J 


(27) 


where  R  =  yj (p2  +  C2  —  l)2  +  4£2  as  before,  r2  =  1  +  (2,  ( tan(0)  =  1,  and 
(P2  +  C2  -  1)  tan(0)  =  2(. 

Sneddon’s  solutions  (equation  5)  can  be  rewritten  more  succinctly  using  the  substitution, 

A  =  and  by  normalizing  the  stress  components  by  the  mean  pressure  prn  =  ppfjzppj,  and  the 
displacement  components  by  the  indentation  depth  e,  resulting  in 


Ut(P,  0  (CJ11-(1~2U)J10) 

e  7r(l  —  u) 

uz(p,  C)  (2(1  —  u)Jq  +  (J®) 

e  7r(l  —  v) 

Trz{p,C)  C^2 

Pm  2 

t'm  ^ 

=  _  J-[2„pJ?  -  C  J\  +  (1  -  2k)  J1]  , 

Pm  " P 

=  “[J?p  -  CM  -  Ji1)  -  (1  -  2k)  Jo1]  .  (28) 

Pm  ^P 

Thus,  the  dimensionless  analytical  solutions  for  displacements  and  stresses  in  an  elastic 
half-space  (equation  28)  are  functions  of  Poisson’s  ratio  v,  and  substantiates  the  choice  of  v  alone 
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to  describe  the  elastic  properties  of  the  half-space  in  the  dimensional  analysis  for  Hertzian  cone 
cracks  developed  in  Gazonas  et  al.  (5);  these  equations  are  corroborated  by  those  found  in 
Barquins  and  Maugis  (27),  with  the  exception  of  the  expression  for  ur,  misprinted  in  both 
Barquins  and  Maugis  (27)  and  Sneddon  (4). 


3.  Peridynamics  Code  Verification 


In  the  previous  section,  we  have  shown  how  to  derive  a  new  explicit  expression  for  J\\  necessary 
for  determining  the  vertical  displacements  in  an  elastic  half-space  (equation  28);  the  solution  for 
the  vertical  displacements, 


MM)  2“s  («-»>  i1  +  Vi  +  MM)))  +  ((<_,Iv/Up,o) 

e  7r 

of  the  half-space  is  illustrated  in  figure  2  (solid  lines),  where  we  verify  a  numerical  peridynamic 
solution  (dashed  lines)  at  dimensionless  depths,  (  =  0,  (  =  1.7,  and  (  =  5.4.  A  100  x  100  node 
grid  spacing  was  used  for  all  axisymmetric  peridynamic  simulations.  The  perfectly  matched 
layer  (PML)  developed  by  Wildman  and  Gazonas  (3)  for  2-D  peridynamics,  is  employed  in  this 
axisymmetric  solution  at  finite  radial  and  vertical  distances  p  =  (  ~  8  within  the  computational 
domain  to  simulate  the  zero  vertical  displacement  condition  that  exists  at  infinity.  The  solution 
for  the  radial  displacements, 


MM)  +  + 

/  -1  \  ?  W'-'/ 

e  7i \y  —  1  )p 

of  the  half-space  is  illustrated  in  figure  3  (solid  lines),  where  we  verify  our  numerical  peridynamic 
solution  (dashed  lines)  at  dimensionless  depths,  (  =  0,  (  =  1.7,  and  (  =  5.4.  The  analytical 
radial  surface  displacements  of  the  half-space  are  generally  non-monotonic  and  the  computed 
peridynamic  results  accurately  reproduce  the  theoretical  trends  as  a  function  of  depth  within  the 
half-space;  we  note,  however,  that  the  computed  radial  displacements  undergo  a  small  spatial 
wrinkle  at  p  =  1  at  the  corner  of  the  rigid  punch,  and  we  surmise  that  this  is  related  to 
singularities  present  in  the  stress  and  strain  fields  that  may  influence  the  non-singular 
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displacement  fields.  Since  we  do  not  observe  such  behavior  at  p  =  1  in  the  computed  vertical 
displacement  field  illustrated  in  figure  2,  it  is  also  possible  that  this  is  a  result  of  a  solution  that 
has  not  fully  converged. 

Using  the  theory  of  the  potential,  Fabrikant  (32;  page  221,  equations  (53)  and  (54))  has 
independently  derived  analytical  solutions  for  both  the  bonded  and  unbonded  punch  problem; 
Fabrikant’s  solutions  are  to  within  machine  precision,  numerically  equivalent  to  the  new  solutions 
uz  in  equation  29  and  ur  in  equation  30.  The  mean  pressure  pm  =  induced  in  the  elastic 

half-space  with  rigid  punch  radius  a  =  5  mm,  indentation  depth  e  =  0.5  mm.  Young’s  modulus, 

E  —  65  GPa,  and  Poisson’s  ratio,  v  =  0.20  is  4.3  GPa. 


Figure  2.  Verification  of  our  peridynamics  code  results  (dashed  lines)  using  the 
analytical  solution  (solid  lines;  equation  29)  for  normalized  vertical 
displacements  ^  in  an  isotropic  elastic  half-space,  v  =  0.2,  due  to  a 
rigid  flat-ended  cylindrical  punch;  the  mean  pressure  induced  in  the 
half-space  is  pm  =  4.3  GPa. 
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Figure  3.  Verification  of  our  peridynamics  code  results  (dashed  lines)  using  the 
analytical  solution  (solid  lines;  equation  30)  for  normalized  radial 
displacements  —  in  an  isotropic  elastic  half-space,  v  =  0.2,  due  to  a 
rigid  flat-ended  cylindrical  punch. 
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4.  Stress  Trajectories 


Stress  trajectories  or  isostatics  (principal  stress  directions  in  the  plane  9  -  constant)  for  the 
problem  of  the  half-space  indented  by  a  rigid  flat-ended  cylindrical  punch  first  appear  in 
Sneddon  (4)  and  can  be  determined  by  solving 


tan(2a)  =  — - — — —  (31) 

(Jr  0~z 

for  a  inclined  at  angles  a  and  a  +  7t/2  to  the  radial  axis,  r.  Use  of  the  auxiliary  relationship  in 
dimensionless  coordinates,  tan(a)  =  d(/dp,  we  arrive  at  a  differential  equation  for  solution  of 
the  isostatics  for  maximum  and  minimum  principal  stress  directions  in  a  half-space  under  the 
action  of  a  rigid  flat-ended  cylindrical  punch: 


d(  ~  crz  (ar  -  a z\ 2 

dp  2 Trz  V  V  2Trz  ) 


(32) 


Substitution  of  stress  components  ar,  az,  and  Trz  from  equation  28  into  equation  32  results  in 


d<  2(pj2°  +  (1  -  2Q Jj  -  (j;  /,  .  (2< p,n  +  (1  -  2 v)4  -  C 

dp  2(  pJl  +V  4(  V(4)2 

To  determine  the  isostatics  for  the  rigid  punch  indentation  problem,  equation  33  (for  the  cr3 
trajectory)  and  its  negative  reciprocal  (for  the  <j\  trajectory)  must  be  numerically  integrated  with  a 
computational  software  program  like  Mathematica  8.0  (35).  The  maximum  shearing  stress 
trajectories  bisect  the  angle  made  between  the  orthogonal  principal  stress  trajectories,  thus,  one 
set  of  maximum  shearing  stress  trajectories  can  be  determined  from  the  equation 


dC  =  q-  1 

dp  q  +  1 

and  an  orthogonal  set  of  maximum  shearing  stress  trajectories  can  be  determined  from 


(34) 
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(35) 


d(  _  q+1 
dp  1  —  q 

where  q  is  given  by  the  right-hand  side  of  equation  33.  The  isostatics  and  maximum  shear  stress 
trajectories  have  been  determined  numerically  using  Mathematica  8.0  (35)  and  these  curves  are 
illustrated  in  figure  4.  In  general,  our  numerically  computed  trajectories  compare  quite  well  with 
Sneddon’s  (4)  hand-drawn  curves,  generated  “by  tracing  two  sets  of  curves  which  are  at  every 
point  tangential  to  the  direction  of  the  two  principal  stresses  in  the  plane.”  Our  <j\  trajectories 
(solid  black  lines  in  figure  4),  however,  tend  to  parallel  the  free  surface  in  regions  outside  the 
contact  area,  a  result  shown  also  in  Barquins  and  Maugis  (27;  page  342,  figure  4),  whereas 
Sneddon’s  curves  (4;  page  38,  figure  6)  are  in  error,  and  orthogonal  to  the  free  surface  in  such 
regions.  Because  our  maximum  shear  stress  trajectories  (dashed  lines)  should  intersect  the  free 
surface,  i.e.  along  the  p-axis  at  45°,  as  they  do  along  the  (-axis,  they  are  in  error  here,  as  we  were 
unable  to  numerically  integrate  equation  35  along  the  p-domain  in  regions  where  the  trajectories 
become  vertical  and  have  undefined  slope,  ^  =  oo;  we  attempted  to  re-couch  equation  35  in 
terms  of  polar  coordinates  to  eliminate  the  numerical  singularity,  but  were  still  unable  to  obtain  a 
solution  in  this  region. 
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p 


Figure  4.  Axisymmetric  isostatics  in  an  elastic  medium,  u  =  0.1679,  indented  by 
a  flat-ended  cylindrical  punch;  a±  stress  trajectories  (solid  black  lines), 
(T3  stress  trajectories  (dashed  black  lines),  and  principal  shear  stress 
trajectories  (orthogonal  dot-dashed  blue  lines). 
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5.  Conclusions 


In  this  report,  we  have  derived  new  expressions  (equation  26)  for  the  ./”'  terms  that  appear  as 
stress  and  displacement  components  in  Sneddon’s  solution  to  the  axisymmetric  indentation 
problem  (1).  The  new  solutions  are  used  to  verify  a  newly  developed  peridynamics  numerical 
code  (2,  3).  The  analytical  solutions  for  the  normal  and  radial  displacements  of  an  elastic 
half-space  under  the  action  of  a  rigid  cylindrical  punch  compare  well  with  the  elastic 
deformations  from  peridynamic  simulations.  Principal  stress  and  shear  trajectories  in  the  elastic 
half-space  under  the  action  of  a  rigid  cylindrical  punch  were  derived  by  numerical  integration  of 
equation  33;  our  results  generally  compare  well  with  Sneddon’s  hand-drawn  curves,  but 
Sneddon’s  oq  trajectories  ( 4 )  are  in  error  near  the  surface  of  the  half-space  and  outside  the  region 
of  indentation. 

Sneddon’s  classic  textbook  (7)  was  published  only  three  years  before  the  publication  of  the 
Bateman  Manuscript  Project  treatise  (25),  so  the  reason  why  Sneddon  was  not  aware  of,  or  did 
not  publish,  any  result  for  the  Jq  term,  given  by  equation  4  (page  101,  equation  (17)  of  that 
treatise),  remains  an  open  question;  the  omission  of  the  Jq  term,  and  hence  expressions  for  the 
vertical  displacements  within  the  half-space,  is  also  apparent  in  more  recent  contact  mechanics 
textbooks  (12,  36).  Finally,  we  derived  a  new  expression  for  the  Jq  term  with  a  relatively 
straightforward  application  of  L’Hospital’s  rule  and  have  shown  that  it  is  analytically  equivalent 
to  the  results  of  Barquins  and  Maugis  (27)  (equation  19),  Gerrard  and  Harrison  (equation  20),  and 
that  found  in  Bateman’s  treatise  (equation  4).  The  new  analytical  results  for  the  remaining  ./”* 
terms  (equation  26)  compare  exactly  to  within  machine  precision12  with  those  from  Sneddon’s 
textbook  (1)  (equation  27)  for  a  number  of  half-space  depths  (  and  radial  distances  p. 


12 


This  result  is  shown  analytically  for  the  Jq  term  in  the  text  and  analytically  for  the  Jq  term  in  the  appendix. 
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Appendix.  Equivalence  of  Sneddon’s  J]  Term  With  Our  New  Derivation 


In  this  appendix,  we  derive  an  analytical  equivalence  between  Sneddon’s  Jq  term  found  on  page 
462,  equation  (56)  of  his  textbook  (f)  with  our  own  derivation,  viz., 


1  — 


(C  -*)  (i  +  tf(p,0) 


(C-i)(H(p,  0-1) 


(A-l) 


Expanding  the  first  term  in  equation  A-l  using  definitions  R  =  y  (p2  +  C2  —  l)2  +  4£2,  and 
(p2  +  C2  —  1)  tan(0)  =  2£  yields 


Jo  =  ^(l-  \/ ^2  +  P2  ~  l)2  +  4C2sin  Q  tan  1  (p2  +  (2  -  1,  ,  (A-2) 

and  the  third  term  in  equation  A-l  using  11(  f).  (J  =  Jl  +  yields 


Expanding  the  term  under  the  square  root  sign  in  equation  A-3  and  taking  the  imaginary  part 
results  in 


J]  =  i  ^1  +  y/ (C 2  +  p2  -  l)2  +  4C2sin  Q  arg  (p2  +  (C  -  f)2)^)  ^  •  (A-4) 

Expanding  the  arg  function  term  in  equation  A-4  results  in 

Jl  =  -p  (i+  \J (C2  +  P2  —  !)2  +  4C2  sin  tan-1  (C2  +  p2  -  1, -2()^  ,  (A-5) 

the  required  equation  A-l  after  factoring  the  minus  sign  out  of  the  arctangent  function. 
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